Elasticity of a system with non-central potentials 
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We derive expressions for determination of the stress and the elastic constants in systems com- 
posed of particles interacting via non-central two-body potentials as thermal averages of products of 
first and second partial derivatives of the interparticle potentials and components of the interparticle 
separation vectors. These results are adapted to hard potentials, when the stress and the elastic 
constants are expressed as thermal averages of the components of normals to contact surfaces be- 
tween the particles and components of vectors separating the centers of the particles. The averages 
require the knowledge of simultaneous contact probabilities of two pairs of particles. We apply the 
expressions to particles for which a contact function can be defined, and demonstrate the feasibility 
of the method by computing the stress and the elastic constants of a two-dimensional system of 
hard ellipses using Monte Carlo simulations. 
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I. INTRODUCTION 



The mechanical response of materials to deformations 
is described by the elasticity theory P|. The simplest 
homogeneous (affine) deformation of a continuum can be 
expressed by a linear dependence of the distorted position 
r on the original position of that point R via relation 



(1) 



where My is a constant tensor 



We will consider both 
two-dimensional (2D) and three-dimensional (3D) sys- 
tems, in which Latin subscripts will indicate Cartesian 
coordinates. (Summation over repeated Latin subscripts 
is implied.) Tensor My can be separated into an identity 
tensor / and a non-trivial part 



(2) 



of the system to assure that the same equation describes 
every internal point, while in the case of inhomogcneous 
system, application of such a deformation to the bound- 
aries assures that the mean deformation is equal to e„ 
i. 

Elastic properties of a condensed matter system de- 
scribe the energetic cost of a deformation. However, real 
system consisting of many moving atoms/molecules can- 
not be simply represented by a strain tensor assigned to 
every point in space. Rather, we can assume that the 
boundaries of such system undergo an affine deformation 
described by Eq. ^ In such a case, the mean free en- 
ergy density, /, which is the free energy F divided by 
the original (unstrained) volume Vo of the system, can 
be expanded in a power series in the strain variables 

/({'7}) = /({O}) + (JtjVij + ^C^jkiVtjVki + ■■■ (4) 



In general, can be separated into a symmetric tensor The coefficients in this expansion are the stress tensor 



representing deformation and an anti-symmetric tensor 
representing rotation. We neglect the rotation and as- 
sume that eij = Cji . While the tensor e.y has a convenient 
meaning in actual experiments, usually the elastic defor- 
mations are formulated in the terms of the Lagrangian 
strain tensor rjij , which defines the change in the dis- 
tances between the points 0: 



M.kMuRkRi = (Ski + 2r,ki)RkRi 



(3) 



In the case of affine deformations the definitions in Eqs. 
are valid for arbitrarily values of eij and rjij . How- 
ever, we will further assume that the deformations are 
small. For a homogeneous continuum, it suffices to ap- 
ply the deformation described by Eq.^at the boundaries 
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aij and the tensor of the (second order) elastic con- 
stants Cijki, characterizing a given material. In the 
case of isotropic pressure p, the stress can be written as 
aij = —p6ij . Elastic constants may serve as an indicator 
of instabilities associated with phase transitions 0, |^ . 

Elastic response of a system to a deformation can be 
determined without actually distorting the system, since 
equilibrium correlation functions contain all the neces- 
sary information. Indeed almost four decades ago Squire, 
Holt and Hoover (SHH) ;6| developed a formalism that 
extended the theory of elasticity of Born and Huang 
to a finite temperature situation and expressed the elas- 
tic properties of a system as thermal averages of various 
derivatives of interparticle potentials. In a certain sense 
the formalism is an extension of the virial theorem Q 
which relates the thermal averages of the products of 
interparticle forces and the interparticle separations to 
the stress tensor. (Similar formalism enables evaluation 
of the elastic properties of 2D membranes in 3D space 
0.) The SHH method is very well adapted for use in 
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numerical simulations in constant volume (and shape) 
ensembles. Other methods, extracting the elastic prop- 
erties from shape and volume changes of systems, have 
also been developed and extensively used 

Usually molecules are not spherically symmetric and 
we may expect interactions that depend on the orien- 
tation of the molecules. The introduction of rotational 
degrees of freedom into the theoretical description of a 
system has an interesting effect on the stress and elastic 
constants. At very low densities (almost ideal gas) the 
rotational degrees of freedom add a large contribution to 
the total kinetic energy of a molecule, but do not con- 
tribute to the pressure. In a thermodynamical treatment 
of stress we are interested in the translational degrees of 
freedom. However, for non-spherically-symmetric poten- 
tials, the particle rotation may have a significant indirect 
contribution even at moderate densities. At larger den- 
sities phase diagram may be strongly influenced by those 
degrees of freedom. The general approach of SHH 6] 
(see also Refs. |^ ^) was applied in a detailed form 
to the case of particles interacting via central two-body 
forces. However, the formalism can be easily extended to 
the systems of particles interacting via non-central po- 
tentials, as will be shown in Section ITU 

Modelling of various systems frequently involves par- 
ticles that interact via hard potentials which are either 
or oo: in fact simulation of 2D hard disk system dates 
back to the origins of the Metropolis Monte Carlo (MC) 
method ■ An obvious reason for the use of such po- 
tentials in simulations is their numerical simplicity. How- 
ever, there are important physical reasons for such mod- 
els: in many situations entropy plays a dominant role 
in physical processes, and the absence of energy scale in 
hard potentials "brings out" the entropic features of the 
behavior. Hard sphere systems have been the subject of 
an intensive research for several decades now (see 13] and 
references therein). They serve as the simplest models for 
real fluids, glasses, and colloids. The phase diagram of 
hard spheres is well known. In 3D this system under- 
goes an entropically driven first-order phase transition 
from liquid to solid phase Q . Elastic constants of such 
solids have been explored in the past 0, 0| . Entropy 
also plays a crucial role in the systems containing long 
polymers, such as gels and rubbers 0,0,^]. Not sur- 
prisingly, hard sphere potentials have been extensively 
used to represent excluded volume interactions between 
the monomers (see and references therein). Kantor 
et al. PH introduced tethering potential, that has no en- 
ergy but limits the distance between bonded monomers, 
to represent covalent bonds in polymeric systems. Such 
hard potential combined with hard sphere repulsion can 
be used to simulate a variety of polymeric systems. Re- 
cently, Farago and Kantor [23| adapted the formalism of 
SHH 6] to hard potentials. This new formalism enabled a 
study of a sequence of entropy-dominated systems, such 
as 2D [23 and 3D |23| gels near sol-gel transition and 
other systems ^2^. 

Orientation of non-spherically-symmetric molecules 



ts a crucial role in the properties of liquid crystals 
For instance, the nematic phase is translationally 
disordered but it has orientational order of the molecules. 
From the early stages of the liquid crystal research it has 
been realized that the entropic part of the free energy 
related to non-spherical shapes of the molecules, by it- 
self, can explain many of the properties of the systems 
j27j |. Not surprisingly, hard potentials were frequently 
used to investigate the properties of liquid crystals. Even 
such simplifications as infinitely thin disks [2^ or in- 
finitely thin rods provide valuable insights into the 
problem. A slightly more realistic picture is provided 
by hard spheroids (3^]. Such simulations were primarily 
motivated by the desire to understand the liquid phases. 
However, two interesting solid phases have been detected: 
both phases are translationally ordered, but only one of 
them has orientational order of spheroids. The orien- 
tational order is absent only when the spheroid resem- 
bles a sphere. Sufhciently oblate or prolate spheroids 
are orientationally ordered in the solid phase. During 
the last twenty years hard spheroids have be studied in 
great detail [3ll. A similar hard potential system that 
is suitable for the study of the liquid crystals is a collec- 
tion of hard spherocylinders (cylinders capped at their 
ends by hemispheres). These molecules have slightly 
more complex phase diagram (which includes smectic- 
A liquid-crystalline phase), and also have been studied 
in great detail ,32]. Like spheroids, they have two solid 
phases. (Spherocylinders do not have a shape resem- 
bling oblate spheroid.) Taken together, spheroids and 
spherocylinders provide a rather coherent picture of in- 
fluence of molecular shape on the phase diagram (see, 
e.g., [13). Hard potentials also have been used in other 
ways to represent non-spherically symmetric molecules 
by combining several spheres or disks into more cornpli- 
cated shapes, such as heptagons 34], or long rods |35l| . 
To make the models more realistic, sometimes attractive 
interaction has been added to the usual hard repulsive 
potential [3^ . 



In Section ^ we present the formalism for soft non- 
central pair potentials. This formalism cannot be di- 
rectly applied to the calculation of elastic constants of 
hard potential systems. In Section IIIII and Appendix 
IXI we show how generalization of the approach used for 
centrally-symmetric hard potentials (23| can be used to 
derive expressions applicable to non-central hard poten- 
tials. In Section IWI we detail the method by which the 
formal results can be applied to hard particles for which 
a "contact function" can be defined. In particular, we 
describe how our formalism can be used for a 2D system 
of hard ellipses. In Section we demonstrate the imple- 
mentation by calculating stress and elastic constants in 
different phases of a system of hard ellipses. 
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II. ELASTIC PROPERTIES FOR SOFT PAIR 
POTENTIALS 

In this section we derive explicit expressions for the 
stress and elastic constant tensors following the method 
of SHH for a more general case. We will consider 
potential energy V which can be expressed as 



(a/3) 



(5) 



where $ is the interaction potential of a pair of 
particles. Greek indices a and /3 denote particles 
(atoms/molecules), and (a/3) denotes a pair of particles. 
The above equation contains summation over all possi- 
ble particle pairs. (In this paper we do not assume sum- 
mation over repeated Greek indices indicating particles.) 
Here r"^ — j./3 _ j.« jg ^^j^g vector connecting two parti- 
cles, while 57" is the orientation of particle a. The two- 
body potential is not necessarily spherically symmetric. 
In fact, we will apply our results to particles that do not 
posses such a symmetry. We denote all two-body po- 
tentials by the same letter $ although nowhere in this 
formal derivation it is required that they should be iden- 
tical for different pairs of particles. (It should be denoted 
$«^(r"^,fl",rj''); however, we omit the superscript of 
for brevity.) From the physical point of view we expect 
that the potential should be rotationally invariant, i.e. 
when the vector r"^ and the orientations of two molecules 
(described by fi" and fl^) perform a "rigid body" rota- 
tion, the interaction energy should not change. In fact 
the symmetry of the stress tensor assumes the presence 
of rotational invariance. However, we do not explicitly 
use this property in the derivation of the following ex- 
pressions. 

Unlike the central force case we will need to use 
both rjij and in the process of derivation. Note that 
in the definition of rjij in Eq. |21 only the symmetric sum 
Vij + Vji appears for i ^ j- Therefore, without loss of 
generality it is assumed that the Lagrangian strain is a 
symmetric tensor (rjij = rjji). From Eqs.[21and|31we find 
that r]ki = ^(cfez -I- e/fc -I- ^ik^u)- For small deformations 
this relation can be inverted to the second order as 



1 



rjki 



(6) 



In statistical-mechanical description of a solid in a 
canonical ensemble we may ask how the free energy F 
of the solid changes when the boundaries of the solid un- 
dergo deformation described by Eq. ^ In calculation of 
such F({r]}) we do not impose any restrictions on the po- 
sitions or orientations of the particles except the change 
in the boundary conditions. The free energy can be ex- 
pressed via the partition function as 



F(M) = -fcTln[ZoZ(M)], 



(7) 



where only the configurational part Z{{ri}) of the par- 
tition depends on the deformation, while the remaining 



("kinetic") part Zq is independent of deformations. We 
note, that in classical physics the details of the inertia 
tensors of the molecules can modify the details of their ac- 
tual motion, but play no role in the statistical- mechanical 
properties of the system. Only the asphericity of the po- 
tential matters. The configurational part 



V({rj}) 

(8) 

where r" and represent the position and orientation 
of particle a and V is the interaction potential, depends 
on the deformation only through the distortion of the 
integration volume V{{i]}) of the possible positions of 
each of the particles. The integration over all possible 
spatial directions of each particle remains unchanged. If 
we formally change the integration variable r" for each 
particle a to the variable R" , which are related by Eq. ^ 
then, the limits of interaction of the new variables will 
correspond to the undistorted volume V^({0}) = Vo, and 
consequently 



„ N N 

z = / n '^^^ / n di^^j'^ 

•J \ 1 ■J \ 1 



-v(Mijij),n\...,Mijflf ,n")/fcT 



A=l " A=l 



(9) 

where J is the Jacobian corresponding to the change of 
coordinates 



J = I det(M)|^ = I det(/ + 27^)1^/2 



(10) 



The deformation now appears as distortion of the coor- 
dinates in the potential V. 

Thus, the stress and the elastic constants can be viewed 
as the first and the second derivatives of the free en- 
ergy density with respect to various rjij . Since in the 
expansion in Eq. 0] always appear pairs of terms such 
as Cii23?7ii'723 + C'ii32?7ii?732 which contain identical rj 
terms, we can choose to define the tensor in a symmet- 
ric form Cijki — Cjiki — Cijik- (An additional symmetry 
Cijki — Ckiij is also evident from the definition of the ten- 
sor.) Strictly speaking, since 7712 = 7721 they should be 
treated as a single variable while taking the derivatives of 
the free energy density. However, terms containing those 
two variables also appear twice in Eq. ^ Thus, one can 
simply treat 7712 and 7/21 as independent variables, and 
symmetrize the results with the interchange of indices 
at the end. Alternatively, one may view each derivative 
9/97712 as ^{d/di]i2 +d/dri2i). Below we always present 
fully symmetrized expressions. 

From Eqs. 01 and [7| we can express the stress tensor 



Vna, 



dF 



kT dZ 

{,,} = {0} ^ 



(11) 



{,,}={o} 
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and the second order elastic constants 



{v}=m 

kT dZ dZ kT d'^Z 



Z"^ dr\r,%n dVij Z drjmndrj^' 



(12) 



in terms of the derivatives of Z . As can be seen from Eq.|^ 
the dependence of Z on the deformation is contained in 
the Jacobian J and in the arguments of the potential V. 
The Jacobian depends directly on riij and its derivatives 
can be easily calculated. In particular, we find (see, e.g., 
Ref. "21) that 



dj 

9mj 



W={o} 



(13) 



9V 



dri,nndr]r 



{n}={0} 



(14) 

Taking the derivatives of V involves the differentiation of 
the potential with respect to My , followed by the differ- 
entiation of Mij with respect to eki using Eq. [3 followed 
by the differentiation of eki with respect to rjmn using 
Eq. El This leads to the following expressions for the 
stress and elastic constants: 



Vna, 



-NkTS. 




(15) 



VoCr 



—y 



R 



,Q/3 



9$ 



Q/3 J 



OR 



R 



:af3 



E 



(75) 



OR 



R2^ 



,7(5 " 



dm" 



— y 



a/3 aD7<5 " 



R^R 



OR 



RZ^Rl^ + 



9$ a$ 



Q/3 r)lS 



a0 ap75 ™ -J' 



dR"^ dR 



R„ R, 



R'^^Rf 




9$ a$ 



dp n7<5 



dRZ^ dRf " ' 



dR^dR'; 



apc,r/ap rn i 



(16) 



ff^Rf, where f"'' is the 



In the above equations we already use the coordinates 
R"'' of the undistorted system to emphasize the fact that 
all the averages are now calculated in the absence of the 
deformation. 

Since {d^/dRf)Rf = 
force between the particles a and /3, we can recognize in 
Eq.^lthe standard virial theorem, although in the text- 
books |8J the derivation is quickly reduced to calculation 
of the (isotropic) pressure p. 

The accuracy of the expressions El and El can be 
verified by reducing the above formulae to the case of 
isotropic central force potential in which case 



/Ri 



(17) 



where prime denotes a derivative of $ with respect to 
inter-particle separation R. Similarly, 



dR,dRj 



, RiRj 
~R2- 



R 



(18) 



After performing these substitutions, we recover the stan- 
dard expressions for central force potentials @. 



III. HARD POTENTIALS 

The expressions for stress and elastic constants that 
have been obtained in the previous section, presumed 
smooth potentials with well defined first and second 
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R 



ap 



a 



is the probability density of particles a and /3 to be at 
particular positions and orientations. The integral in 
Eq. ^1 simply represents the average of n^^R"^ over 
all possible contacts between two particles weighted with 
proper probability densities of those contacts. (Since the 
probability density changes from a finite value, when the 
particles are almost in contact, to zero, when the par- 
ticles overlap, we need to consider the case when s"^ 
approaches from the positive side.) Thus, if we denote 



, a/3 



0+), we can find the hard potential 



limit by replacement 



FIG. 1: Two "hard" particles at a close approach. R"*^ 
is a vector connecting their centers, and s"'^ is the minimal 
separation of the surfaces. 



dK 



Q/3 



-kTA. 



a/3 



(21) 



derivatives. There is a certain difficulty in translating 
the expressions obtained for soft potentials to a hard po- 
tential situation. For instance, the term {d^/dR"^)R"'^ 
in Eq. E| looks poorly defined for a hard potential that 
changes on contact between and oo. However, we note 
that the derivative of the potential really originates from 
the derivative a[e-*(^"'''""'"'')/'=^]/ai?f The latter, 
is a derivative of a step function that changes between 
and 1, when the potential changes between oo and 0. 
This observation, has been used in Refs. 0, to de- 
rive simple expressions for the pressure of hard sphere 
system. A detailed and rigorous description of the var- 
ious aspects (including calculation of pressure) of inter- 
action of hard convex solids can be found in Ref. [s^. 
A typical interaction between two hard particles (e.g., 
ellipses in 2D or ellipsoids in 3D) is represented by Fig. 
n which depicts a close approach of two such particles, 
when R"'' is the vector connecting their centers, while 
s"^ denotes the minimal distance between them. Gradi- 
ent of e-*(^°^'"°'"'')/'='^ with respect to R"^, when 
and fl^ are kept fixed, is simply n°'^S{s°'^) where n"'^ is 
unit vector perpendicular to the surfaces at the point of 
contact, pointing from a to (3. Thus we can make the 
substitution 



Rf\^-kT dR"(iR'3(if7°rff7'^ 



\dR^^^^' 
<''<5(s"'')i?"'^P(R", r!", R^, nl^) 



(19) 



where 



P(R",f7",R'5,17'^) 

1 r ^ 



^<°f> (20) 



By substituting the result into Eq. E| we obtain the fol- 
lowing expression for the stress in the case of hard po- 
tentials: 



kT 



= -N6,. 



\E 



AfRf 



(22) 



(a/3) 



The substitution appearing in Eq. 1211 cannot be used 
to calculate the elastic constants in Eq. ^1 since the lat- 
ter expression contains a double summation (over pairs 
(a/3) and (7(5)) of the derivatives of the potentials. Such 
summation includes a term where (a/3) = (7*5). A direct 
substitution of Eq. [2] would lead to the appearance of 
the term [5{s°'^)]^ which causes the expression to diverge. 
(This is a true divergence, rather than some mathemati- 
cal subtlety of approaching the limit of hard potentials.) 
However, Eq. ^| also contains the second derivative of 
<i>, which becomes poorly defined in the hard potential 
limit. We shall show in the Appendix 1X1 that the sum of 
the apparently divergent and poorly defined terms has a 
well defined hard potential limit. In fact this sum can 
be transformed into an expression in Eq. IA5I which in- 
cludes only first derivatives of potentials and products 
of the first derivatives of potentials of different particle 
pairs. While the resulting expression in Eq. IA5I looks 
more complicated, it can be simply transformed (using 
Eq. I21|) into an expression for the elastic constants of 
hard particles: 



A = l 



6 



kf 



(a/3) (75) 
(a/3) <7«> 

+ ^ E E (^'' + ^"') ^"''^ + (^"'' + ^T) RZ^'Rf 



(a/3) 

Af i?^'^ + A^'^ijf ) <5„ + (A^:f i?f + Af % + (Af i?f + Af i?f ) 5,„„}. (23) 



Each of the terms in the above equation corresponds to 
some particular case of contact between pairs of parti- 
cles. E.g., term of the type {A"'^ R"^) in the last sum 
corresponds to a contact between a pair of particles a 
and /3, and consequently the contribution of that sum 
is proportional to N. The sum with the prefactor i 
contains terms corresponding to three touching particles: 
E.g., the term {/S.<^ A'f R'^'^ R'^'^) corresponds to the sit- 
uation when particle /3 touches particles a and 7 simul- 
taneously. The number of such contacts at any given 
moment is also proportional to TV. The second and the 
third lines of the equation have different averages 
corresponding to the number of pairs of contacts that 
might appear in the system. However, we note that in 
both sums taken together we always have differences of 
the type {Af Rf){A'^^Rl') - {Af Rf A^J R^') . This 
term vanishes if the contact between the pair {a(3) is un- 
correlated with the contact between the pair {jS) . This 
happens when the two pairs are outside the correlation 
distance. Consequently, only pairs of contacts close to 
each other will contribute, and therefore the total contri- 
bution of those terms is proportional to TV. 



IV. APPLICATION TO HARD ELLIPSE 
SYSTEM 

Equations [221 and [221 enable calculation of the stress 
and elastic constants of a system consisting of hard par- 
ticles, provided we are able to calculate thermal averages 
of the type {A"^ R"^) . This expression depends on the 
probability density of particles being in contact. For the 
calculation of stress we need only the probability den- 
sity of a single contact between a pair of particles, while 
the calculation of elastic constants involves the proba- 
bility density of two such contacts happening simultane- 



ously. It is natural to use MC method "i^, 'il', "i^l to 

evaluate averages of this type. MC simulation of hard 
potentials is particularly simple since no energy scale is 
present, and every elementary move of a particle is ei- 
ther accepted without a need to calculate the Boltzmann 
factor, or results in a forbidden configuration, and is, 
consequently, rejected. Application, of these procedures 
requires a method to identify intersection of two hard 
particles. Such methods have been found both for 2D 
ellipses [S^ and for 3D ellipsoids 0,||4^]- In Appendix[Bl 
we explain in detail the case of ellipses which is used in 
the current work. Here, we will consider a slightly more 
general case when a simple function can be defined that 
identifies the contact between two particles. 

Consider a function 4" that depends on the positions 
and orientations of two particles, that vanishes when the 
particles touch each other and is positive when the parti- 
cles are separated. (The formalism can be trivially gen- 
eralized to the case when the the function has some other 
non- vanishing constant value at the contact.) Consider 
a case when the orientations of both particles are fixed, 
and we explore positions where the function vanishes. 
Figs.[2t and[2|3 depict 2D cases of one fixed ellipse, while 
another (identical) ellipse is rotated by a fixed angle and 
is shown in a variety positions where they touch. The ra- 
tio between the major and minor semi-axes of the ellipses, 
a and b, respectively, are different in both pictures. The 
thick line traces the possible positions of the center of the 
moving ellipse. Note that the shape of such a contact line 
depends on the degree of elongation of the ellipses and 
on their relative orientation. Those are the positions that 
are relevant for the calculation of the average {A"^ R"^) . 
In 2D this is a line, while in 3D this is a surface. 

The direction of the force, normal to the contact plane, 
is also normal to this surface. Thus, assuming that 
^' is a sufficiently smooth function, we can calculate 



(a) 




(b) 



FIG. 2: (a) Center of slightly elongated ellipse {a/b = 3) 
tilted by 45° with respect to an identical (vertical) ellipse cir- 
cumscribes an "oval" trajectory, depicted by the thick line, 
when the contact point moves along the ellipse, (b) Cen- 
ter of strongly elongated ellipse {a/b = 14) rotated by 90° 
with respect to an identical (vertical) ellipse circumscribes a 
"rounded square" trajectory, depicted by the thick line, when 
the contact point moves along the ellipse. 

n"^ = (9^'"'3/9i?f'')/|V«'"''|, where the gradient (and 
the partial derivative) with respect to R"^ is taken when 
the orientations of the particles are fixed. In thermal av- 
erage we need to calculate 

{AfRf) = J dn'^dnP X 

J dS°'^nfRfP(R",n",Ti^,n'^), (24) 

which involves the integration along the contact 
surface (or line) S°'^ of the probability density 
P(R",l]",R'3,fi/'), defined in Eq. EOI of two particles 
being in those positions and orientations. During a MC 
simulation such an event strictly never occurs. We can 
replace the integration along the surface, by an integra- 
tion inside a thin shell of thickness t along the contact 
surface. In such a case dS^^P « {dV/t)P = dp/t, where 
dV is the volume element and dp is the probability of the 
center of a particle being within a shell, at some partic- 
ular area. Note, that the thickness of the shell does not 



FIG. 3: Line of * = and * = 1000 for a pair of elhpses 
with aspect ratio E = a/b — 2 with their axes rotated by 45° 
for a function defined in Appendix|B| Note that the thickness 
of the area between two lines of fixed "l/ varies slightly. 

have to be constant, but can vary from place to place on 
the contact surface. In fact, we can define the shell as 
corresponding to all positions for which < 'i'"^ < 5*0, 
where ^'o is some fixed number. Fig. |31 depicts such a 
shell corresponding to two values of '^"^ , for a function 
defined in Appendix ^ If is small enough, we can 
determine the local thickness of the shell from the rela- 
tion « iV^'^'^li. Substituting the values of t and n"'^ 
expressed via function into Eq. and noting that 
approximate expressions mentioned in this paragraph be- 
come exact for vanishing ^Pg, we arrive at the expression 

(Afpf) = hm / dQ'^dn'' / ^^Rfdp 

S(*o) 

(25) 

(AfRf) = lim -L / ^Rf\ (26) 

In these equation S'('I'o) in the integral and in the thermal 
average denote a shell defined by the limit ^E'o on the 
function VP"''. 

In the numerical calculation of the stress, the limit 
in Eq. |^ is not easy to implement: Using a large ^'o 
leads to an inaccurate answer, while using a small ^'o 
leads to a small number of "almost contact" events, and, 
consequently, to large statistical errors. One may try 
considering a numerical extrapolation to ^'o = by mea- 
suring the stress for a sequence of decreasing ^qS. How- 
ever, the events for smaller values of ^'o are also con- 
tained in set of the events for larger ^fgs. It is diffi- 
cult to extrapolate such correlated sets of data points. 
The independence of the data points can be achieved by 
calculating a sequence of values of the stress for "con- 
tact shells" defined by '^ap located in a sequence of seg- 
ments [0, *o), [^-0, 2*o), • • • [K^o, {K + l)*o), • • • (Kis 
integer). Values of (7^ now can be conveniently extrapo- 
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lated to their "real" values. A similar method has been 
used by Farago and Kantor to calculate the stress 
and elastic constants of hard sphere solids. 

The terms in the expressions for elastic constants in- 
cluding two pairs of particles can be similarly handled. 
One simply uses ^'oi and 4*02 to define two shells, S'(\E'oi) 
and £'(^'02) 7 respectively, each corresponding to a partic- 
ular contact, and considers the events which occur when 
both pairs of particles are within their respective shells 
simultaneously. E.g., 



\ m "'n 2 "'7 / 



lim 



1 / d^'^^ 0d^^' s\ 



OR 



(27) 



S(*0l),S(*02) 



Compared with the case of the stress, the numerical eval- 
uation of the limit where the thickness of the shells van- 
ishes presents an even bigger numerical problem, since 
the probability of two contact events is very small. Nev- 
ertheless, this can be handled in a similar way, by con- 
sidering a 2D array of segments of the type {[K'^q, [K + 
l)^'o), [L^'o, (i + l)*o)} {K and L integers) and obtain- 
ing the values of the various parts of the elastic constants 
by extrapolating the 2D surface to its "real" value of van- 
ishing contact layer thickness. 



V. RESULTS OF SIMULATIONS 

We used the method developed in this work to calcu- 
late the elastic properties of 2D hard ellipse system as 
a part of a study of its phase diagram [43 • Here, we 
briefly demonstrate the usefulness of the method. As 
in any hard particle system, temperature plays no role, 
since the interactions have no "energy scale." The tem- 
perature appears only as a multiplicative prefactor in the 
free energy F and in Eq. |^for the stress and Eg. 1231 for 
the elastic constants. The results depend on the density 
of the particles and their size and shape: we character- 
ized the system by the number of particles per unit area p 
and by the sizes of the major and minor semi-axes, a and 
6, of the ellipses. Frequently, reduced density p* = Apab 
is used. The maximal possible (close packed) p* is inde- 
pendent of the aspect ratio E = a/b oi the ellipses and is 
equal to 2/-\/3 « 1.155. It should be noted [s^ that for 
every fixed a and b there is an infinity of possible (equally 
dense) close packed states which are obtained by orient- 
ing all ellipses in the same direction and packing them 
into a (distorted triangular) periodic structure. 

The system of hard disks {E ^ 1) has been extensively 
studied. For p* > 0.91 it forms a periodic 2D solid — a 
triangular lattice. The correlation function of atom po- 
sitions of such a solid decays to zero as a power law of 
the separation between the atoms Such behavior is 
usually denoted as quasi-long-range order. At the same 
time the orientations of the "bonds" (imaginary lines con- 
necting neighboring atoms) have a long range correlation 



[?n| . The system is liquid for p* < 0.89, i.e it has no long 
range order of any kind. At the intermediate densities 
the system is probably hexatic 0| — a phase with alge- 
braically decaying bond-orientational order, but without 
positional order. (However, even very large scale simu- 
lations 49] have difficulties in distinguishing the hexatic 
phase from coexisting solid and liquid phases.) From 
the point of view of elasticity theory, all three phases 
are isotropic, i.e. their second order elastic constants are 
determined by two independent constants. Aspect ratio 
E 1 oi the ellipses adds an additional order param- 
eter — their orientation. E.g., for £' = 4 a system of 
ellipses forms isotropic liquid for densities p* ^ 0.8. For 
larger densities the ellipses in the liquid become oriented 
("nematic phase"). Finally, at p* « 1.0 the system be- 
comes a solid of orientationally ordered ellipses |50i| . For 
weakly elongated ellipses, we expect the particles to re- 
main orientationally disordered in the entire liquid phase, 
and with increasing density to go (possibly via hexatic 
phase) to a crystalline state in which the ellipses remain 
disordered. The 3D analog of such a state is called plastic 
crystal 51]. (Presence of such a phase in almost circu- 
lar ellipses {E = 1.01) was observed by Vieillard-Baron 
[3^.1 With increasing density an additional phase transi- 
tion will bring the ellipses into an orientationally ordered 
state. Phase diagram which includes such a transition 
between two solid phases for 3D hard ellipsoids has been 
studied by Frenkel et al. 30]. 

As a test of our formalism we studied a case of mod- 
erately elongated hard ellipses with E = 1.5. We consid- 
ered system consisting ol N k, 1000 ellipses contained in 
a 2D rectangular box whose dimensions were chosen as 
close as possible to a square. Periodic boundary condi- 
tions were used. The ellipses were initially placed on a 
distorted triangular lattice, commensurate with the di- 
mensions of a closed packed configuration corresponding 
to this particular aspect ratio E and the particular ori- 
entation of the ellipses. In this Section we describe only 
the cases when the initial orientation of the major axes of 
the ellipses was taken to be perpendicular to one of the 
axes of symmetry of the ordered crystal drawn through 
neighboring particles. We first performed an equilibra- 
tion run at constant pressure, in order to allow the system 
to reach an equilibrated state with respect to the orien- 
tational, as well as the translational ordering. The ori- 
entational order parameter, the box dimensions and the 
density were monitored during this equilibration run. A 
MC time unit in the equilibration run consisted of + 1 
elementary moves, one of which, on the average, was a 
volume change attempt, and the rest were particle move 
attempts. A particle move attempt involved choosing a 
particle randomly, and attempting to displace and rotate 
it simultaneously by an amount chosen from a uniform 
distribution. The move was accepted if the displaced par- 
ticle did not overlap in its new position and orientation 
with other particles. The volume change attempt was 
identical to that described in jij. The width of the dis- 
tributions corresponding to the particle moves and the 




p p 



FIG. 4: Pressure p (open circles connected by solid line) and 
the elastic constant C1212 (inverted open triangles connected 
by dashed line) of hard ellipse system with aspect ratio E — 
1.5 in the units of kT/Aab as functions of the reduced density 
p*. The error bars of the pressure are significantly smaller 
than the symbols denoting the data points. 



FIG. 5: Elastic constants Cim (open circles connected by 
solid line), C2222 (asterisks connected by dotted line) and 
C1122 (full triangles connected by dot-dashed line) of hard el- 
lipse system with aspect ratio E = 1.5 in the units of kT/Aab 
as functions of the reduced density p* . The error bars of all 
the data points are slightly smaller than the symbols denoting 
them. 



parameter of the volume change were chosen so that the 
average success rates of both types of MC moves were 
about 50%. 

The use of the constant pressure simulations at the 
equilibration stage enabled us to determine the equilib- 
rium box shape for isotropic stress tensor (pressure) con- 
ditions. At most densities the final configuration had 
orientationally disordered ellipses, and we could easily 
verify that the final configurations were independent of 
the specific choice of the starting configuration. How- 
ever, at extremely high densities, approaching the close 
packed density, even after long equilibration the state re- 
sembled the starting state of orientationally ordered el- 
lipses. Typically, several millions of MC time units were 
required to reach equilibrium for a given pressure. Upon 
completion of equilibration, we switched to constant vol- 
ume simulations, during which the stresses and the elas- 
tic constants were evaluated. Very long simulation times 
(about 10^ MC time units) were required for accurate de- 
termination of these constants, because their calculation 
depends on extremely rare events of two pairs of parti- 
cles simultaneously touching each other. The range of 
"contact shells" was chosen in such a way that even in 
the most remote shell the separation between the parti- 
cles was significantly smaller than their mean separation. 
The statistical accuracy of the elastic constants was eval- 
uated by comparing the results of independent runs. 

In a 2D system which has a reflection symmetry with 
respect to either x or y axis, the elastic constants with 
an index appearing an odd number of times (such as 
C1112) must vanish. Indeed, our simulations showed that 
these quantities vanish within the error bars of the mea- 
surement. The system still may have four unrelated 
elastic constants Cim, C2222, (^1122 and Cyivi- For a 



system with quadratic symmetry the number of inde- 
pendent elastic constants reduces to three. Such sys- 
tems are frequently characterized by their bulk modu- 
lus and two shear moduli, \x\ = C1212 — P and /i2 = 
^(C'liii — (^1122) — P- For an isotropic system jii = ^2 
and, therefore, there are only two independent constants. 
(A system with six-fold symmetry is isotropic as far as 
the elastic constants are concerned.) Figs. ^and|Sldepict 
the pressure and four elastic constants for several values 
of the reduced density p* . One can see that for densities 
p* < 1.086, Ciiii and C2222 practically coincide, indicat- 
ing that these systems are at least quadratic. Further- 
more, the identity Cnn — C1122 = 2C1212, i-e. Hi = /i2, 
is found to hold within a few percent for these values 
of p* . Consequently, for these densities the system is 
isotropic from the point of view of the elastic proper- 
ties. Figs. Et, b and c depict the system in that range 
of densities: Figs.lHt- and b represent states with vanish- 
ing shear moduli, and neither of them exhibits transla- 
tional order of the ellipse centers. However, while Fig.^t 
represents a state with all the characteristics of a liq- 
uid, the state in Fig.O) is characterized by slowly decay- 
ing bond orientational order, possibly indicating hexatic 
phase. At a slightly higher density. Fig. ^ represents 
a plastic solid: while the ellipses are randomly oriented, 
the system exhibits long range bond orientational order, 
and algebraically decaying positional order of the parti- 
cles. The particles occupy, on the average, positions of an 
undistorted triangular lattice although some undulations 
are apparent. This is characterized by two coinciding 
positive shear moduli. 

When we approach within few percent the close packed 
configuration, corresponding to the two largest densities 
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in Figs.0]and|Slthe prolonged relaxation process does not 
change the preferred orientation of the ellipses and the 
distortion of the lattice. It would be reasonable to con- 
clude, that at such high density we finally arrived at the 
orientationally ordered elastic solid. The isotropic elastic 
symmetry no longer holds. We checked and found that 
almost all stability criteria ^ , indicating the sign of 
the energy change upon small deformation, are positive. 
However, one of the shear moduli, namely /zi, is slightly 
negative, although within one standard statistical devia- 
tion from zero. This, may either indicate that we are in 
an unstable state, or that there is a continuum of equi- 
librium states with different mean orientations of ellipses 
and, correspondingly, different dimensions of elementary 
cell. 



VI. DISCUSSION 

We extended the formalism of SHH to the case of 
systems interacting via non-central two-particle poten- 
tials. In its form represented by Eqs. El and El the 
formalism can be used to study properties of molecular 
systems. This is particularly true for highly non-spherical 
organic molecules, and various soft condensed matter sys- 
tems. The adaptation of the expressions to hard poten- 
tials (Eqs. 1221 and 123 produced slightly more compli- 
cated expressions. However, the simplicity of hard poten- 
tial systems provides excellent insights into the entropy- 
dominated systems. We demonstrated the usefulness of 
the formalism by presentin g; so me results of our study 
of the hard ellipse system '45|. Measurement of sev- 
eral order parameters, and the correlation functions is 
not always sufficient to determine the nature of phases. 
For systems of moderate size, it maybe even difficult to 
distinguish liquid from a solid. Measurement of elastic 
constants provides an additional, very important tool for 
assessing the nature of the state of the system. In partic- 
ular, the elastic constants may indicate the presence of 
instability, even when prolonged equilibration does not 



change an existing state. Following the indications of in- 
stability at high densities, we are currently performing 
extensive study of equilibrium states at these densities. 

While we worked on a 2D example, the method can 
be equally well applied in 3D for such systems as hard 
ellipsoids or spherocylinders. 
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APPENDIX A: REGULARIZATION OF 
DIVERGING AND POORLY DEFINED TERMS 

Section IIIII outlines the procedure for transition from 
expressions for soft potentials to expressions for hard po- 
tentials. The procedure relies on the fact, that despite the 
change in the potential between and oo, the expression 
for stress and some terms in the expression for for elas- 
tic constants contain a derivative of Boltzmann factor. 
Since the latter changes between 1 and 0, its derivatives 
involve (5- function of the distance between the particles, 
leading to a finite thermal average. Even the terms con- 
taining contacts between two distinct pairs of particles, 
and, consequently, two i5-functions of different variables, 
produce a finite results. However, the third line in Eg. 1161 
includes a term containing a product of two derivatives of 
the same pair of particles. A substitution of Eg. 1191 into 
such an expression would lead to appearance of squared 
(5-function and divergence of the entire term. Eg. 1 161 also 
contains second derivative of which becomes poorly 
defined in the hard potential limit. Here we show that 
sum the two "problematic" terms described above has 
a well defined hard potential limit. Let us consider one 
particle pair (a/S): 



1 a$ 9$ 



a/3 r>a/3 



R^R 



kTdR'^idRf dRZ^dRf j " ' 

Qi g-$(R°'',n°,n'')//cT 



dR"df7° dRf^dnf^— — -s RfRfp(R'',n'',R'^,n'^) (ai) 



In the internal integral in Eg. lAll we can replace d/dRf by d/dR^ , since the exponent only depends on R'^ — R". 
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FIG. 6: Typical equilibrium configurations of slightly eccentric {E = 1.5) hard ellipse system at several densities. Only part of 
the system is shown. All the pictures show the same (partial) volume of the system; they differ only in the reduced density p* : 
(a) orientationally and translationally disordered liquid at p* = 0.881; (b) liquid with a high degree of bond-orientational order 
at p* — 1.029; (c) plastic solid with long-range bond-orientational order, and quasi-long-range translational order consisting of 
rotationally disordered ellipses at p* — 1.086; (d) solid of orientationally ordered ellipses at p* — 1.117. 



Following that, we perform integration by parts in which boundary term vanishes and arrive at 



kT / dR'^dn" / dR^dn^— ^ -^[RfRfp(R",n",R^,n^)] 

J J dR'if dR^ ' 

= - I dR-^dn" f dR^dn^-^ g-*(R°^o^a'')/fcT ( s,^^Rfp + s,,RfP + RfRf^] . 



(A2) 



The probability density P in the above expressions was defined in Eq. 1201 Variable appears in every potential that 
depends on some R''^, and therefore, the derivative of P with respect to R^ in the last term of the above expression 
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can be expressed in the following form: 



dP 1 
Z 



dR 



N 



A = l 



dm 



e 5*(°'3> 



A = l 



,7/3 



e ^{c'fi) 



(A3) 



Consequently, the term in Eg. lAll becomes: 



9$ 



'■j 



Si?; 



-E 



9^ d'^ r^cPr^o^P 

.dRZ?dRf ' , 



(A4) 



(This answer is slightly non-symmetric — this reflects the fact that we manipulated the integral over the variable R'^, 
rather that R". In the latter case instead of partial derivative with respect to Rj^ we would have obtained a partial 
derivative with respect to R'^ ■ ) When this result is substituted in the expression for the elastic constants in Eg. 1161 
we get 



VoC, 



ijmn 



AkT 



-y 

(a/3) 

E E 



iR 



dRf ' 



dRf 



R 



afi 
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(75) 
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a/3 0715 



4fcr 



(a,3) {1^) 



8kT ^ ^ 



R^, R, 



■R°''^R- 
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dR^ 



R 



dRZ^ dR. 



RZ^R''^ 
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a/3 p7'5 
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9$ 
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dR 
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-5, 



a/3 
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9$ 
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9$ 



9i?: 



a/3"''J 



a/3 



dR 



,a/3 "-J 



7 1 "-m "-^ 



(A5) 



Since the above expressions contain only the first derivatives of the potentials (or products of such terms for non- 
identical pairs of particles) we can use Eas.ll9landl21lto calculate the elastic constants for hard potentials. 



APPENDIX B: OVERLAP FUNCTION OF TWO whether two such objects overlap. Vieillard-Baron de- 
ELLIPSES rived a contact function for the 2D case of ellipses fs^l- 

His function enables a reasonably fast determination of 

A crucial difficulty in the simulations of hard ellipses 
in 2D and hard spheroids in 3D is the determination 



13 



the presence of an overlap. It can also be used to de- 
termine the direction of the inter-ellipse force and, con- 
sequently, can be used in determination of stresses and 
elastic constants. In this Appendix we describe this func- 
tion and its properties. 

Consider two identical ellipses, a and (3, whose major 
and minor semi-axes are a and 6, respectively. Let their 
centers be separated by vector R"^, and one of them 
be rotated compared to the other one by the angle 0. 
Projections of R"^ on the directions of major and mi- 
nor semi-axes of ellipse a will be denoted and i?^2 ' 
respectively. Then we can define a contact function 



4{hl - 2,hp){hl - ZK) - (9 - Khpf, (Bl) 



where for each ellipse we define 

K = l + G-[{Rl1/aY + {Rl^/b)\ (B2) 



and 



G 



I sm 

a , 



(B3) 



Function depends on the distance between the el- 
lipses and on their orientation. The number of overlap 
points of two ellipses can vary between and 4. The 



necessary and sufficient condition for ellipses to have no 
intersection points is > and at least one of the 
two functions ha and /i^ be negative. The proof of this 
statement can be found in Ref. [s^l- At the contact 
yj,a/3 _ ^]^g region where the ellipses do not inter- 

sect the function ^'"'^ has no extremum points; thus if 
vl'"'^ is sufficiently small, i.e. < ^'"^ < ^'o we can be 
assured that two ellipses are close to each other. 

When the ellipses degenerate into a circles, i.e. a = 
b, the contact function becomes only a function of the 
distance between the centers of the circles i?"'^ 



^"/3 ^ 3 



-4 



(B4) 



Note, that the gradient of this function at the distance 
of contact between the circles has a rather large value 
768/a, so even if we choose — 8 the circles will be 
about 1% of the radius away from each other. 

The shell between the lines ^""^ = and 5'"'^ = has 
a uniform thickness for circles, but its thickness varies 
with position in ellipses as can be seen in the Fig. 13 
When the aspect ratio of the ellipse becomes much larger 
than unity the thickness becomes strongly variable which 
will adversely influence the accuracy of Monte Carlo sim- 
ulations. 
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